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A numerical method to calculate optical conductivity based on a pump-probe setup is presented. 

Its validity and limits are tested and demonstrated via concrete numerical simulations on the half- 
filled one-dimensional extended Hubbard model both in and out of equilibrium. By employing 
either a steplike or a Gaussian-like probing vector potential, it is found that in nonequilibrium, 
the method in the narrow-probe-pulse limit can be identified with variant types of linear response 
theory, which, in equilibrium, produce identical results. The observation reveals the underlying 
probe-pulse dependence of the optical conductivity calculations in nonequilibrium, which may have 
applications in the theoretical analysis of ultrafast spectroscopy measurements. 
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I. INTRODUCTION 

Time-resolved spectroscopy has been employed exten¬ 
sively to investigate the dynamic properties of materials 
with strong correlations, e.g., (quasi-)one-dimensional 
(ID) Mott insulators such as bis(ethylendithyo)- 
tetrathiafulvalene-difluorotetracyanoquinodimethane 
ET-F2TCNQip— and halogen-bridged transition-metal 
compounds — as well as correlated systems with higher 
space dimensions, such as cuprates^^S and mangan- 
ite o^^d^ . As a promising time-resolved spectroscopy 
technique, ultrafast pump-probe optical measurement is 
able to unravel the complicated entanglement of various 
degrees of freedom in a correlated system to some extent 
by resolving their contributions separately at multiple 
time scalesi^. For experiments on electronic systems, 
time-resolved optical conductivity is an essential quan¬ 
tity, from which the dynamic properties of charges and 
other related degrees of freedom can be analyzed; see 
the discussions in Refs. [H andH 

Generally, in order to calculate a dynamic response 
function, the corresponding correlation functions at dif¬ 
ferent times should be obtained. In nonequilibrium, the 
task becomes more demanding since the system is evolv¬ 
ing and the time-translation invariance is absent. Con¬ 
sequently, the correlation functions with two time vari¬ 
ables cannot be reduced to functions with only one time 
variable, i.e., the temporal distance in the case of equi¬ 
librium. Theoretically, there are various ways in micro¬ 
scopic model calculations to describe the dynamic re¬ 
sponse of correlated systems in nonequilibrium. One pos¬ 
sible way is based on the nonequilibrium Green’s func¬ 
tions: by using the Kadanoff-Baym or Keldysh formal¬ 
ism, methods that have been applied to correlated sys¬ 
tems in equilibrium, such as dynamic mean-field theory, 
can be adapted to nonequilibrium on a formal level (see, 
e.g.. Ref. [la and references therein). Another method 
begins with wave functions; this method is used pri¬ 
marily within closed systems. The time-dependent wave 


functions can be secured using various numerical meth¬ 
ods, such as exact diagonalization and the density-matrix 
renormalization-group method. Thereafter, the time- 
dependent expectation values of observables and vari¬ 
ous temporal correlations can be readily obtained, which 
enables further calculations on the out-of-equilibrium 
dynamic response. Specifically, with regard to time- 
dependent optical conductivity, several related but dif¬ 
ferent schemes have been employed in similar models, in¬ 
cluding the nonequilibrium generalizations of the Kubo 
formulapi^, linear-response theoryiiil^, or a direct calcula¬ 
tion on dynamic current-current correlations for transient 
statesi^. However, the underlying characters and validity 
of these approaches have not yet been fully addressed as 
far as we know. This is also a practical problem for the 
purpose of analyzing data from (sub-)THz spectroscopy 
measurements, where the temporal resolution has been 
pushed up to the order of a few tens of femtoseconds^. 

In this paper, we suggest that the way in which tem¬ 
poral correlations are incorporated can be crucial in 
the study of nonequilibrium dynamics. Inspired by the 
pump-probe setup in the experiments, we propose a 
method to calculate the optical conductivity. Its va¬ 
lidity and limits are demonstrated via detailed numer¬ 
ical simulations on the ID half-filled extended Hubbard 
model, a prototype of strongly correlated electronic sys¬ 
tems. A comparison between our method and the ex¬ 
isting ones shows that by adopting two different limits 
on the probe pulse, our method can reproduce the re¬ 
sults of two types of linear-response theory. We note 
that in some nonequilibrium situations, the two sets of 
results on the time-dependent optical conductivity can 
be quite different. Theoretical analysis combined with 
numerical simulations elucidates the difference and con¬ 
nections between these various methods to a satisfactory 
level. They also raise the issue of the probe-pulse depen¬ 
dence in nonequilibrium analysis, which has been ignored 
in previous studies. 

The rest of the paper is organized as follows. In 
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Sec. m after the model description and a critical re¬ 
view of the relevant existing methods, our approach is 
presented. The following section provides a theoretical 
analysis on the validity of the approach, with numerical 
demonstrations concerning the options for probe-pulse 
parameters. Section ITVl contains comprehensive compar¬ 
isons and analysis. The conclusion is drawn in Sec. [V] 


II. MODEL AND METHODS 

The Hamiltonian we use in the discussion is the well- 
known ID extended Hubbard model at half-filling: 

H = —tfi 'y ^ H.c.^ 

i,<7 

+ (na-0 

+ vY^{ni-l){n^+l-l), ( 1 ) 

i 

where c\ ^ (ci.a-) is the creation (annihilation) operator 

of an electron with spin a at site i, nt^a = cl^Ci^a, 
n-i = + th is the hopping constant, and U and V 

are on-site and nearest-neighbor Coulomb repulsion in¬ 
teractions, respectively. In the following, we set th and 
as energy and time units and take the Plank constant 
h = 1, and the elementary charge e = 1. 

We restrict our discussions to zero temperature. Before 
presenting our method in detail, we first present some 
existing methods that are relevant to our discussions. 

In the linear-response theory, the induced current due 
to the application of an external perturbation field E{t) 
is given by (in the ID case) 

{Jit'))= f a{t\t)E{t)dt, (2) 

J —CO 

where the response function <y{t', t), which is probe-field- 
independent, measures the current response at time t' 
with respect to the perturbation of a unit pulse at time 
After the response function is determined, the in¬ 
duced current for any given perturbation field E{t) in 
the linear-response regime can be calculated accordingly 
by Eq. (g]). 

The situation can be simplified in equilibrium, where 
due to the presence of time-translation invariance (with¬ 
out consideration of the probe field), the two-time re¬ 
sponse function a{t', t) can be reduced to a single-variable 
function as a(t' — t). Then the Fourier transformation on 
Eq. (21) simply produces 

j{uj) = a{uj)E{uj). (3) 

(j{uj) represents the optical conductivity, which measures 
the current response of the system in frequency space. 
It is worth noting that although the conductivity criui) 
can be evaluated or measured by jiuj)/E{ui) as indicated 


in Eq. m (see additional discussions in Sec. dmi), it 
reflects the property of the system in equilibrium and 
does not depend on the details of the probe field. It is well 
known that cr(w) can be calculated microscopically by 
the Kubo formula (KF)^ in terms of the current-current 
correlations. The regular part (excluding t he p ossible 
singularity at w = 0) is given by (e.g., in Ref. 1^ : 

c^reg(w) = e*“'*(0|[i(s), j(0)]|0) ds = ^x(w), 

(4) 

where |0) is the ground state of the system, L is the 
lattice size, and the electric susceptibility x(w) is defined 
as 

/ -(-oo 

e“'*(-f)6»(s)(0|[j(s),j(0)]|0)ds. (5) 

-OO 

The current j{s) in Eqs. ([3]) and ([5]) is the Heisenberg 
representation of the current operator j, which reads 

j = -ithYyA<y^i+i,^ - H-C.]. (6) 


In nonequilibrium, on the other hand, since time- 
translation invariance does not hold in general, a full two- 
time response function is required in order to describe the 
linear response of the electric current with respect to a 
time-dependent probing field Eit). However, for another 
related and important quantity, i.e., the time-dependent 
optical conductivity that we will mainly focus on, there 
is no unique definition if the time-translation invariance 
does not hold in the response functio n^^i^^i^^ . Following 
Ref.M throughout the article we define the conductivity 
as: 

a{uj,t) = [ (T(t-I-s, ds, (7) 

Jo 

where we have explicitly introduced a spectral broad¬ 
ening factor rj and the time width tm for the Fourier 
transformations that will be used in the numerical sim¬ 
ulations. In practice, tm is the maximum of the probe 
duration time, which in THz spectroscopy is usually sev¬ 
eral hundred times larger than the microscopic time unit 
(e.g., tl'^ here). 

From Eq. it is clear that if we intend to calculate 
a{ui,t) by investigating the current response to a weak 
probe field analogous to Eq. which will be developed 
as the pump-probe method later in the article, a highly 
narrow probe pulse is preferred. However, different from 
the equilibrium case, we find that the response character 
of the system in nonequilibrium depends on the form of 
the probe pulse that we choose, which will be detailed 
in Sec. lYl Here, we first review the relevant theoretical 
formulas related to a{uj^t). 

In the linear-response regime, rigorous analysis pro¬ 
duces the response function a{t',t) as (with restriction 
t’ > t)^ 


a{t',t) = 


{^(t')\rm'))+X{t\ndt'‘ 


, ( 8 ) 
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where the two-time susceptibility is 

x{t\n = ( 9 ) 

and in the diamagnetic term, r = th'^(c\+i rjCi,a + 

Z,(7 

H.C.), which is the stress tensor operator. The inter¬ 
action representation of an operator {t') is defined as 
U\t',t) B where is the time-evolution 

operator in the absence of probing perturbations^. Con¬ 
sistent with this, the time-dependent wave function |'0(t)) 
follows the nonequilibrium time evolution of the (closed) 
system without the probe field. It is clear from the ex¬ 
pression that in order to obtain the value of a{t',t), the 
current-current correlations (commutators) between the 
intermediate time and the ending time t' have to be taken 
into account. We refer to this as nonequilibrium linear 
response (NLR). In the next section, we will show that 
the Fourier transformation of a{t', t) according to the def¬ 
inition (O reflects the current response with respect to 
the unit electric pulse at time t. 

In the literature, we notice that another kind of gener¬ 
alized Kubo formula in nonequilibriuin (NKF) has been 
used^ii^: 


ReCTreg(wji) =—plm 

(jJL/ 


ds 




x{m\[j'{t+s),j\t)]\m)]- ( 10 ) 


A careful analysis shows that the NKF cannot be ob¬ 
tained from the response function in Eq. ([5]) (with the 
diamagnetic term dropped) by using the definition of 
Eq. ([7|). It seems that another ’’response function”, de¬ 
noted as cr{t',t) [t' > t) here, is used: 


a{t',t) = - 


Xit'\t)dt'‘ 


( 11 ) 


where similarly 

x{t'\t) = -i 0 {t"-t)m)\[j'{n,j\mm)- ( 12 ) 


We refer to this as a variant of NLR (VNLR). Eirst, it is 
easy to recognize that the NKE in Eq. m is just the real 
regular part of the Fourier transformation of the second 
term in Eq. (HU with respect to t' — t. Furthermore, the 
difference between NLR and VNLR actually comes from 
the second term in the response functions (Eqs. ([5|) and 
(HU), where under the integration over the intermediate 
time t" € (t, t'), the two-time commutators of the current 
operator take place at t" to t', and t to t", respectively. 
Note that t and t' are the starting time and ending time 
of cr{t',t) it' > t), respectively. That is to say, different 
from NLR, d'{t',t) takes into account the current-current 
correlations (commutators) between the starting time t 
and the intermediate time. The above description can be 
read from Fig. [TJ In the next section, we will show that 
the Fourier transformation of d{t' ,t) according to the def¬ 
inition © reflects the current response with respect to a 
delta vector potential applied at time t. 


NLR 



t t" ... ... t' 


VNLR 



t t" ... ... t' 


FIG. 1. (Color online) A schematic diagram of the current- 
current correlations considered in NLR and VNLR, respec¬ 
tively. 


We argue that the NLR method, which has been de¬ 
rived rigorously from linear-response theoryi^, contains 
more complete information about the temporal correla¬ 
tions than VNLR (or equivalently NKF) with respect to 
the definition of aiujfi) in Eq. The observation can 
be easily seen in Fig. [TJ reminding us that aiujfi) is the 
Fourier transformation of (j{t',t) with respect to t' — t 
(with t fixed). The fact suggests that VNLR can be re¬ 
garded as a kind of partial summation on the temporal 
correlations of NLR. It should be emphasized that for the 
equilibrium case, the two produce identical ct(w) due to 
the time-translation invariance. Then an issue could be 
raised regarding the validity and limitation of VNLR: in 
what situation and to what extend can VNLR remain a 
good approximation? Does it impose any specific restric¬ 
tion on the probing fields? To address this issue, we now 
present our method. 

Inspired by the experimental pump-probe setup and 
the theoretical works on the coherent dynamics in the 
BCS model^i^, we propose a method to calculate the 
optical conductivity of the Hubbard model both in and 
out of equilibrium as follows. 

First, one notes that an external spatially homoge¬ 
neous electric field applied along the chain in the Hamil¬ 
tonian © can be incorporated via the Peierls substitu¬ 
tion in the hopping terms as^: 

^ H^H{t). (13) 

With the knowledge of the wave function ip{t) under the 
action of A(t), the temporal evolution of the expectation 
value of the current operator {j{t)), defined as {j{t)) = 
('*/’(^)|j(^)l'0(O)i can be readily obtained. Note that with 
the existence of A(t), the current operator j becomes 
time dependent: 

jit) = = -tthY^[e^^^''^cl^c,+xa - H.c.]. (14) 

In equilibrium, we set the vector potential A(t) to be 
a weak probe pulse denoted as Api.obe(i), and the current 
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induced by it as (jprobe(t))- Here the time-dependent 
wave function describes the time evolution of the 

system under the influence of Hpi.obe(t) starting from the 
ground state. Identical to Eq. the optical conductiv¬ 
ity can then be calculated by 

\ _ jpioheji^)/L _ _ Jprobe(^) _ 

E'probe(^) d” H^)Z/Hpi.obe(^) 

where jprobe{w) and Hprobe(w) are the Fourier transfor¬ 
mations of (jprobe(t)) and Hprobe(t), respectively. Note 
that numerically, a damping factor is also introduced 
when the Fourier transformations are performed, as in¬ 
dicated in Eq. 0. The same rj in the denominator of 
Eq. (fT^ is necessary to distinguish the Drude component 
in the spectral weight at w —>■ 0, which will be discussed 
later in more detail in Sec. El 

It is worth to stress here that Eq. m can be regarded 
as a practical definition of the conductivity in equilib¬ 
rium: < 7 ( 0 ;) obtained in this way does not depend on the 
details of the probing held we choose, which has been 
verihed in the numerical simulations. We are now in the 
position to generalize the algorithm to nonequilibrium, 
where the system can be driven either by a pump pulse 
or a quench. The key is that in order to identify the 
response of the system with respect to the later probe 
pulse, a subtraction is necessary, i.e., two successive steps 
are involved (each step follows a similar process to that 
described above) in order to calculate the optical con¬ 
ductivity in nonequilibrium^l. First, a time-evolution 
process that describes the nonequilibrium development 
of the system in the absence of a probe pulse is evalu¬ 
ated, and we obtain {j(t)). Secondly, in the presence of 
an additional probe pulse, we get (jtotai(t))- The sub¬ 
traction of (j(t)) from (jtotai(t)) produces the required 
(jprobe(t)), i-e., the variation of the current expectations 
due to the presence of the probe pulse. Then analogous 
to Eq. m, the optical conductivity in nonequilibrium is 
proposed to be 


/ , \ _ Jprobe(^; tprobe) 

p,obej = ,(^ + i^)2,ylp,obe(w)’ 


(16) 


where tprobe is the probing time. This pump-probe-based 
method is abbreviated as PP. The main reason for con¬ 
sidering the PP method seriously is that for us it seems 
to be a reasonable numerical approach to the ultrafast 
spectroscopy experiments. 

A few final words are in order regarding the numerical 
method: to trace the temporal evolution of the system, 
we employ the time-dependent Lanczos method to eval¬ 
uate |'0(t)). The key formula is 


M 


where e; and \4>i) are eigenvalues and eigenvectors of the 
tridiagonal matrix generated in the Lanczos iteration, re¬ 
spectively, M is the dimension of the Lanczos basis, and 
6t is the minimum time stem More details of the algo¬ 
rithm can be found in Ref. l28l. 


III. THE INFLUENCE OF THE PROBE PULSES 
- A THEORETICAL ARGUMENT 


In the PP method, analogous to experiments, it is im¬ 
portant that the strength of the probe pulse should be 
weak enough not to disturb the system’s dynamics qual¬ 
itatively. Meanwhile, the influence of the forms of the 
probe pulse is also crucial in nonequilibrium. This point 
will be made clear in the following arguments. 

We begin with the specific form of Eq. ([2]) in ID: 


(jprobe(t')) = / cr(t', t")^^probe(t") dt", (18) 
Jto 


where we have assumed that the probe electric held is 
turned on after to with the restriction of t' > to- The 
Fourier transformation on (jprobe(tO) gives 


Jprobe(^) 



I'thi 

= / (jprobe(t'))e*"‘' dt' 

Jto 

a(t',t")-Eprobe(t")e'‘"‘' dt"dt' 


plM r^M 

= / fT(t',t")e“^* ■* ^dt' 

Jto Ut" 

ptM 

= a(c^,t")7?probe(t")e“* dt", 

Jto 


Eprobe(t")e*-‘'' 

(19) 


where Im is the maximum evolution time and ^ tg. 
We can see that in order to obtain (T(w,t) [assuming t S 
(to, tAf), and in this section we use t to denote the prob¬ 
ing time], the most convenient choice of Eprobe in theory 
is a (5 pulse with the form of ifpi.obe(t") ~ S{t" — t), and 
the desired cr(a;,t) is simply proportional to jprobe(<j-’). 
The strategy is different from the usual treatment for 
equilibrium in the textbooks, e.g., in Ref. where a 
monochromatic perturbation held E{t) ^ jg 
from t —>■ —oo. In the temporal gauge, the favorable 6 
electric held corresponds to a step-like vector potential, 
i.e., Aprobe(t") = Ao,step where 0(t") is the Heav¬ 

iside step function (the effect of ramping, i.e., the gradual 
increase of Aprobe from zero to a hnite value, instead of 
a jump, will be discussed later in this section). We call 
this kind of PP way with a step-like vector potential a 
PP-step. From Eq. (HSl) it is easy to convince oneself 
that under the step potential, Eq. (HU) exactly produces 
a(uj,t) in NLR. This means that by applying step-like 
vector potentials at time t, we can extract the exact value 
of by simply looking at the time evolution of the 

induced current (jprobe(t))- Thus the connection between 
PP-step and NLR can be established. This proposition 
will be verihed by later numerical simulations. 

However, in the ultrafast spectroscopy measurements, 
THz ac probe pulses, which usually come from the same 
source of pump pulses, are employed instead of electric 
d-like helds. It seems that a more realistic approach for 







5 


(b) ,! 

h 

- 

f ^probe^^^ 

ik. 

r 

w 
;! ’ 

i[ * -PP-Gaussian with 

} - - - PP-Gaussian with 2*t 



(d) 

1 

1 Eprobe(^) 

A 

1 

-PP-ramp with td 

U - PP-ramp with 2*td 

\\ -PP-step 



FIG. 2. (Color online) Schematic illustrations of the vector 
potentials and the corresponding electric fields for the Gaus¬ 
sian pulse [(a) and (b)] and the steplike pulse [(c) and (d)] 
with various widths. 



FIG. 3. (Color online) The calculated Re ( 7 ( 0 ;) in equilib¬ 
rium (zero temperature) for the Hamiltonian ([T]) with differ¬ 
ent fd,probe . Parameters: L = 10, [/ = 10, and V = 4.5. (a) 
By PP-Gaussian with oiprobe = 10, and Aq, probe = 1.0 x 10“®. 
(b) By PP-ramp with Ao,step = 1.0 x lO”'^. Note that the 
difference is negligible between the results of td,probe = 0.2 
(red curves) and those of td,probe = 0.02 (blue curves). 


a probe pulse could 


.^probe(^ ) —^0,probe exp (t 
X COS [Wprobe (<" “ t)] 


tf /2I 


2 

d,probe 


( 20 ) 


where a Gaussian-like envelope around t is used with 
id,probe to characterize the temporal width of the probe 
pulse, and Wprobe is the central frequency. Following the 
previous conventions, we call the PP method with the 
Gaussian-like vector potential PP-Gaussian. An illustra¬ 
tion of Gaussian pulses with different time widths td is 
shown in Fig. [2][a). 

Let us examine what approximation has to be made 
in order to justify Eq. (ITSl) in the PP-Gaussian. Sup¬ 
pose that the pulse is activated during [t,t + At], i.e., 
Eproheit") 7 ^ 0 only when t" e [t,t + At]. Gonsequently, 


(jprobe(^O) = 0 when t' < t. From Eq. (IT^ we have 

~ (7(cJ, t)£/probe(^); (^1) 

which is nothing but Eq. ([TC|). The key approximation 
lies in the last step of Eq. (ED), where a{uj,t") is ap¬ 
proximated by <7{uj,t). It requires a narrow probe pulse, 
i.e., td,probe should be small enough for nonequilibrium in 
order to obtain It is noted that even in equilib¬ 

rium, the width of the probe pulse can also influence the 
resulting a{uj), as shown in Fig.O In Fig.[3][a), the zero- 
temperature optical conductivity cr(u;) in equilibrium for 
the Hamiltonian (IT|) on the L = 10 lattice is calculated by 
PP-Gaussian under different values of td,probe- It shows 
that when td,probe < 0.2, the results are converged. The 
reason is as follows: in PP-Gaussian, the incoming pho¬ 
ton frequency is broadened into a Gaussian-like distribu¬ 
tion, with the variance of 1/t^ around Wprobe; in order to 
cover a large enough w-regime, a sufficiently small td,probe 
is preferred. This observation explains the large devia¬ 
tion for td,probe = 0.6 at low frequencies when Wprobe is 
set to be 10 in Fig. [31[a). 

To make the information complete, we also investi¬ 
gate the ramping effect in PP-step, where a schematic 
illustration of Api.obe(i) with a finite increase in width 
(called the PP-ramp) can be found in Fig. [21[b) (simi¬ 
lar to a half-cycle pulse in Ref. IH). For the Hamilto¬ 
nian with identical parameters used in the previous 
PP-Gaussian calculation, we see that the results converge 
as td,probe < 0.4, as shown in Fig. EKb). The noticeable 
deviation for td,probe = 0.6 at high frequencies is due to 
the exponential decay of Api.obe(<j-') in the frequency space 
with the decay constant proportional to p,.oi 3 e- 

The above discussion tends to suggest that with small 
enough td,probe, the results of PP-Gaussian and PP-step 
would be the same. We have shown that this is actually 
the case in equilibrium. However, there is a critical dif¬ 
ference between them. Mathematically, note that in the 
PP-Gaussian, unlike the PP-step, Aprobe(t") —0 as t" —>• 
±oo, and consequently, the temporal integration over the 
electric field gives /// E{t") dt" « /+“ E{t") dt" = 0, 
which might be true in experimental situations^. The 
induced current (jprobe(tO), shown in Eq. (flSl) . is pro¬ 
duced in PP-Gaussian by a temporally localized pertur¬ 
bation in the range of [t, t + At] (outside the range, both 
the electric field and the vector potential are practically 
zero). It is therefore reasonable to suggest that when 
the Gaussian-like pulse is used, only the correlations be¬ 
tween the variant t' and (roughly) fixed t (probe time) 
are taken into account in the calculation of cr(w, t), which 
coincides with the feature of VNLR when the temporal 
correlations are concerned (see Fig. [T|). 

It is instructive to consider a limiting case in which 
we simply put Wprobe = 0, and td,probe -t 0, i.e., using a 
(5-like vector potential. If we put an electric field corre¬ 
sponding to the S vector potential (similar to a monocy- 
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FIG. 4. (Color online) Re(T(ct)) in eqnilibrinm (zero tempera¬ 
ture) for the Hamiltonian (IT|) obtained by the KF, NLR, and 
PP methods in (a) the open boundary condition (open BC) 
and (b) the periodic BC. Parameters: L = 10, U = 10, and 
V = 4.5. 

cle pulse in Ref. into Eq. ([2]) by using the response 
function in NLR, and if we perform a partial integration, 
we find that the resulting (jprobe(^O) is proportional to 
x(<^,t) (here we just focus on the current-current corre¬ 
lation term). Employing its Eourier transformation into 
Eq. (1161) , what we finally obtain is merely the expression 
of the NKE in Eq. (fTT|l . The loop between PP-Gaussian 
and VNLR thus can be closed. From the above analysis, 
we propose that the PP-step and PP-Gaussian can be 
identified with the NLR and VNLR, respectively. This 
proposition will be examined in the next section numer¬ 
ically. 

In the following simulations, we set tm 200, M = 30, 
St = 0.02, and rj = 1/L, with the lattice size L = 10. PP- 
Gaussian and PP-step (as the limiting case of zero-width 
PP-ramp) are employed, compared with the results of 
NLR and VNLR. In the PP-step, Ao.step = 1-0 x 10“^; 
in the PP-Gaussian, Aq, probe 10“®, td,probe = 0.02, 
and Wprobe = 10. Only the real parts of the optical con¬ 
ductivity are displayed. Despite double temporal evolu¬ 
tions being involved, the performance of the PP methods 
can match those of (V)NKF in speed, and they are less 
memory-demanding because no integration is requiredi^. 


IV. RESULTS AND COMPARISON 

In this section, working on the ID half-filled extended 
Hubbard model ([T]), we compare the numerical results of 
the optical conductivity both in and out of equilibrium 
between various methods. 

We first discuss the equilibrium case. Figure |4] shows 
the results in equilibrium (zero temperature) under open 
and periodic boundary conditions (BCs), obtained by the 
KF, NLR, and PP. There are no differences between the 
results of the PP-step and PP-Gaussian. The NKE, as 
a generalization of the KF to nonequilibrium, produces 
data identical to the KF in equilibrium. We can see that 
the results obtained by the three methods coincide well 
in the open BC, while in the periodic BC, deviations 
occur in the vicinity of a; = 0. As shown in Fig. HKb), an 
identical peak around zero frequency is produced both 
by the PP and NLR, while it is absent in the KF result. 


The small peak is known to originate from the nonzero 
charge stiffness (Drude weight) in the periodic BC for 
finite systems, even when the system is in an insulating 
phase. The reason is as follows. The charge stiffness D 
is related to the singularity of Recr(a;) at w = 0 as22i2^: 

ReCT(w) = 2'itD5{lij) -\- ReCTreg(w), (22) 

and here 

D = (l/2L)[(r)+Rex(c^ = 0)], (23) 

where x(w) is defined in Eq. ([S]). We can see that the 
Drude weight component comes from the net contribu¬ 
tion of the diamagnetic term (r) and the current-current 
correlations at w = 0. For an insulator, D might still re¬ 
main finite especially in the periodic BC when the system 
size is small, though it should vanish in the thermody¬ 
namic limit. It has been confirmed by the sum rule check 
(not shown here) that for our model at half-filling with 
L = \Q, D has a sizable value in the periodic BC, while 
in the open BC the value is negligibly small. The in¬ 
troduction of the rj factor in the numerical calculations 
expands the singularity at a; = 0 into a small peak, as 
shown in Fig. SKb). Here we can see that the contribu¬ 
tion to the conductivity from the Drude weight can be 
captured both by the NLR and PP methods, while the 
KF method, which only takes into account the regular 
part of the current-current correlations in Eq. (|3]), fails 
to produce this feature. 

We are now ready to proceed to the nonequilibrium 
case. As a demonstration, three well-attended situations 
are considered here. The driving force to nonequilibrium 
in the first case is a pump pulse with a form similar to 
Eq. (1201) (formally just replace ’’probe” with ’’pump”), 
and we call this situation the Pump case. In the second 
case, we apply a temporary (5-like electric pulse, which 
can be simulated by introducing a phase shift in 

the hopping terms of Eq. where 9{t) is the Heav¬ 

iside function. The strength of the pulse is controlled 
by (j). In our calculations we choose 4> = tt, an extreme 
case that is known as 7r-quench. The third case is the 
[/-quench, where a sudden change of U takes place at 
t = 0. This situation can be realized, for example, in 
cold atoms^. 

Figure [S] shows the results of Re(T(u;,t) at / = 10 for 
these three setups with various methods. First note that 
the NKF method only includes the regular part of the 
current-current correlation term (see Eq. (HOI) ). As a re¬ 
sult, Re(T(a;,t) equals to zero at w = 0 exactly. Away 
from the vicinity of w = 0, the results of the NKF 
and VNLR are consistent. This is not surprising since 
the NKF is just about the regular part of VNLR. More 
importantly, from Fig. [Sj we can see that in all three 
cases, including pump and quench, the results between 
PP-Gaussian and VNLR (left column) and PP-step and 
NLR (right column) coincide well with each other. Thus 
the proposition in the previous section of the equivalence 
among these groups is verified numerically. We make 
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FIG. 5. (Color online) The time-dependent optical conduc¬ 
tivity Recr(aj,t) at t = 10 for the Hamiltonian ([T]) obtained 
using different methods in the periodic BC, where t is the time 
delay between the probe and pump (quench) time. In (a) and 
(d), a pump pulse is applied with parameters Ao.pump = 0.2, 
uipump = 6.29 (resonant frequency), and fa,pump = 0.5. In (b) 
and (e), a 7r-quench is applied. In (c) and (f), a [/-quench 
with the U value changed from 10 to 4 is applied. The blue 
solid lines represent the results of the VNLR in the left col¬ 
umn ((a), (b), and (c)) and NLR in the right column ((d), (e), 
and (f)), respectively, while the red dashed lines represent the 
results of the PP-Gaussian and PP-step in the two columns. 
The insets in the left column display the details in the vicinity 
of oj = 0 obtained by NKF, VNLR and PP-Gaussian. Model 
parameters: L = 10, [/ = 10, V = 3. 


it clear that the NKF and NLR, as being widely used 
in the dynamic response study, actually correspond to 
different limits of the probe in the PP method, which 
either has a Gaussian-like vector potential with a vanish¬ 
ing width, or a steplike one. The former with a tiny 
width as small as fd, probe = 0.02 actually can be re¬ 
garded as an approximation of the temporal (5-like po¬ 
tential, which has been numerically confirmed in our cal¬ 
culations (not shown here). Though the electric fields 
are both highly temporally confined in these two cases, 
one essential difference is the persistence (or temporal 
locality in the opposite sense) of the vector potential 
in the subsequent time evolutions, which results in the 
(non)vanishing of the time integral of the electric field 
over the probing time. A possible implication of the 
above observations is that the time-dependent conduc¬ 
tivity or other similar transport quantities in nonequi¬ 
librium obtained in experimental measurements, might 
depend on the details of the probe if the features of the 



FIG. 6. (Golor online) The time-dependent optical conductiv¬ 
ity Re(T((u,t) for the Hamiltonian ([T]) obtained by NLR and 
VNRL in a [/-quench case, where U is switched from 0 to 2 
at [ = 0. The probing time is set at t = 0 in (a) and t = 5 in 
(b), respectively. Model parameters: L = 10 and V = 0. 


measuring procedures are essentially described by the PP 
method. Theoretical investigations should take into ac¬ 
count the relevant experimental setup sufficiently. In the 
present THz spectroscopy measurements, the duration 
time of probe pulses is usually hundreds of times larger 
than the microscopic time unit, which may vindicate the 
approximation of PP-step or NLR, where the character¬ 
istic duration time of the probe is much longer than the 
time scale of the dynamics in which we are interestedi^. 
Conversely, i.e., if the probe time scale is much smaller 
than the latter, VNLR or PP-Gaussian might be a good 
starting point for analysis. 

It may be worthwhile to make some supplementary 
comments on the physics in Fig. [5] before closing the sec¬ 
tion, though it is not the main subject of the present 
study. In the Pump case, a prominent Drude-like struc¬ 
ture of the transient conductivity in the low-frequency 
regime can be observed, associated with the suppression 
of the main absorption peak around w = 6.3. It is a 
typical phenomenon due to photoinduced carriers. More 
discussions on this issue in terms of this model can be 
found, e.g., in Ref. [l^- The 7r-quench corresponds to 
a population inversion in the kinetic energy^, where a 
temporally enhancement of charge mobility is also ob¬ 
served. The U -quench is one of the examples of interac¬ 
tion quench^^. In our case, with U being switched from 
10 to 4, part of the interaction energy is released into the 
system with the consequent enhancement of the kinetic 
energy, which also results in substantial changes in the 
optical spectrum. 

Additionally, as shown in Fig. [SJ a similar develop¬ 
ment of a(u!,t) between the left and right columns can 
still be recognized regardless of the notable difference in 
quenches. This is consistent with the previous statement 
in Sec. El that the VNLR is a kind of partial summa¬ 
tion on the full temporal correlations in NLR. In the 
Pump case, the results of VNLR and NLR are closer 
to each other— compared with the cases of quenches, 
which might be due to the fact that the Hamiltonian only 
changes temporarily with the application of the pump 
pulse, after which it goes back to its earlier form; in 
both TT-quench and [/-quench, the Hamiltonian under¬ 
goes a substantial change after the quench time. Nu- 



















































merical checks in the periodic BC up to L = 14 and the 
open BC in L = 10 show similar results. Thus we tend 
to conclude that the deviation between VNLR and NLR, 
which is due to different temporal correlation counting, 
persists with the increase of the system size. Further, we 
note that the deviation of VNLR from NLR can be quite 
pronounced if the probing time is set to be close to the 
pump or quench moment before the relaxation fully com¬ 
mences. This is because, in contrast with NLR, VNLR 
only considers the temporal correlations starting from the 
probing time t (see Fig.[T|). Figure [5] shows RetT(w,t) for 
the ID half-filled Hubbard model (V = 0 for the Hamil¬ 
tonian ([T])) in a typical U-quench case, where the interac¬ 
tion U is switched from 0 to 2 at t = 0. The probing time 
is set to be t = 0 and t = 5 in (a) and (b), respectively. 
Note that at 17 = 0, the ground state of the system is in 
a metallic state, and with any finite U a charge gap is 
opened and the system becomes insulating. In Fig.lHKa), 
we can observe that the resulting ReCT(w, 0) from VNLR 
shows a strong reminiscence of the initial metallic state 
with a prominent Drude-like peak appearing near zero 
frequency, while the result from NLR, which takes into 
account more information about the subsequent nonequi¬ 
librium evolution, not only shows the depression of this 
peak, but it also produces additional structures in the 
higher-frequency regime. However, when the probing 
time is set at t = 5, by which time the relaxation be¬ 
tween the kinetic energy and the interaction energy is 
largely finished, the results from VNLR and NLR bear 
more qualitative similarities. 


V. CONCLUSION 

In this paper, we introduce a numerical method to cal¬ 
culate the time-dependent optical conductivity based on 
a simulation of the pump-probe experimental setup. Ac¬ 
cording to the property of the probe pulse, two different 
approaches are discussed, i.e., PP-step and PP-Gaussian, 
using either a step-like or a Gaussian-like probing vector 
potential. We find that the two approaches can be identi¬ 
fied with nonequilibrium linear-response theory and one 
of its variants, respectively. Thus the connections be¬ 
tween these various methods and their differences are 
clarified. Numerical verification is given by employing 
the method systematically in the ID half-filled extended 
Hubbard model, both in and out of equilibrium. The 
latter includes three well-attended cases: the pump, tt- 
quench, and [/-quench. In the numerical results, the 
probe-pulse dependence in nonequilibrium is especially 
significant in the quenches. We suggest that the nature 
of the probe pulses should be taken into account in the 
analysis of ultrafast THz spectroscopy measurements. 

ACKNOWLEDGMENTS 

The authors thank Janez Bonca, Denis Golez and Dirk 
Manske for helpful discussions. H.L. acknowledges sup¬ 
port from the National Natural Science Foundation of 
China (NSFC) (Grants Nos. II3254I7, 11474136, and 
11174115), and the open project of State Key Labora¬ 
tory of Theoretical Physics (SKLTP) in ITP, CAS. T.T. 
acknowledges support by the Grant-in-Aid for Scien¬ 
tific Research (26287079) from MEXT and the Strategic 
Programs for Innovative Research (SPIRE) (hpl40215, 
hp 150211), the Computational Materials Science Initia¬ 
tive (CMSI). 


* luht@lzu.edu.cn 

^ H. Okamoto, H. Matsuzaki, T. Wakabayashi, Y. Taka- 
hashi, and T. Hasegawa, Phys. Rev. Lett. 98, 037401 
(2007). 

^ S. Wall, D. Brida, S. R. Clark, H. P. Ehrke, D. Jaksch, 
A. Ardavan, S. Bonora, H. Uemura, Y. Takahashi, 
T. Hasegawa, et al., Nature Phys. 7 , 114 (2011). 

^ M. Mitrano, G. Cotugno, S. R. Clark, R. Singla, S. Kaiser, 

J. Stabler, R. Beyer, M. Dressel, L. Baldassarre, D. Nico- 
letti, A. Perucchi, T. Hasegawa, H. Okamoto, D. Jaksch, 
and A. Cavalleri, Phys. Rev. Lett. 112 , 117801 (2014). 

S. Iwai, M. Ono, A. Maeda, H. Matsuzaki, H. Kishida, 

H. Okamoto, and Y. Tokura, Phys. Rev. Lett. 91, 057401 
(2003). 

® H. Matsuzaki, M. Iwata, T. Miyamoto, T. Terashige, 

K. Iwano, S. Takaishi, M. Takamura, S. Kumagai, M. Ya- 
mashita, R. Takahashi, Y. Wakabayashi, and H. Okamoto, 
Phys. Rev. Lett. 113 , 096403 (2014). 

® D. Fausti, R. Tobey, N. Dean, S. Kaiser, A. Dienst, 
M. Hoffmann, S. Pyon, T. Takayama, H. Takagi, and 
A. Cavalleri, Science 331 , 189 (2011). 


^ D. Nicoletti, E. Casandruc, Y. Laplace, V. Khanna, C. R. 
Hunt, S. Kaiser, S. S. Dhesi, G. D. Gu, J. P. Hill, and 
A. Cavalleri, Phys. Rev. B 90, 100503 (2014). 

® W. Hu, S. Kaiser, D. Nicoletti, C. Hunt, I. Gierz, M. Hoff¬ 
mann, M. Le Tacon, T. Loew, B. Keimer, and A. Cavalleri, 
Nature Mater. 13, 705 (2014). 

® S. Dal Conte, L. Vidmar, D. Golez, M. Mierzejewski, 
G. Soavi, S. Peli, F. Banff, G. Ferrini, R. Comin, B. M. 
Ludbrook, L. Chauviere, N. D. Zhigadlo, H. Eisaki, 
M. Greven, S. Lupi, A. Damascelli, D. Brida, M. Capone, 
J. Bonca, G. Cerullo, and C. Giannetti, Nature Phys. 11, 
421 (2015). 

M. Rini, N. Dean, J. Itatani, Y. Tomioka, Y. Tokura, R. W. 
Schoenlein, A. Cavalleri, et al., Nature 449, 72 (2007). 

T. Li, A. Patz, L. Mouchliadis, J. Yan, T. A. Lograsso, 

I. E. Perakis, and J. Wang, Nature 496, 69 (2013). 

J. Orenstein, Phys. Today 65, 44 (2012). 

J. Orenstein and J. S. Dodge, Phys. Rev. B 92, 134507 
(2015). 

D. Nicoletti, M. Mitrano, A. Cantaluppi, 
and A. Cavalleri, ArXiv e-prints (2015), 

arXiv: 1506.07846 [cond-mat.supr-con] 



9 


H. Aoki, N. Tsuji, M. Eckstein, M. Kollar, T. Oka, and 
P. Werner, Rev. Mod. Pliys. 86, 779 (2014). 

G. De Filippis, V. Cataudella, E. A. Nowadnick, T. P. 
Devereaux, A. S. Mishchenko, and N. Nagaosa, Phys. Rev. 
Lett. 109 , 176402 (2012). 

D. Rossini, R. Fazio, V. Giovannetti, and A. Silva, Euro- 
phys. Lett. 107 , 30002 (2014). 

Z. Lenarcic, D. Golez, J. Bonca, and P. Prelovsek, Phys. 
Rev. B 89, 125123 (2014). 

R. Fukaya, Y. Okimoto, M. Kunitomo, K. Onda, 
T. Ishikawa, S. Koshihara, H. Hashimoto, S. Ishihara, 
A. Isayama, H. Yui, and T. Sasagawa, Nature commun. 
6, 8519 (2015). 

R. Kubo, J. Phys. Soc. Jpn. 12, 570 (1957). 

G. D. Mahan, Many-particle physics (Springer Science & 
Business Media, New York, 2013). 

M. Eckstein, M. Kollar, and P. Werner, Phys. Rev. B 81, 
115131 (2010). 

M. Eckstein and P. Werner, Phys. Rev. Lett. 110, 126401 
(2013). 

T. Papenkort, V. M. Axt, and T. Kuhn, Phys. Rev. B 76 , 
224522 (2007). 


A. P. Schnyder, D. Manske, and A. Avella, Phys. Rev. B 

84, 214513 (2011). 

D. Golez, J. Bonca, L. Vidmar, and S. A. Trugman, Phys. 
Rev. Lett. 109 , 236402 (2012). 

H. Lu, C. Shao, J. Bonca, D. Manske, and T. Tohyama, 
Phys. Rev. B 91 , 245117 (2015). 

P. Prelovsek and J. Bonca, in Strongly Correlated Systems, 
Springer Series in Solid-State Sciences, Vol. 176, edited by 
A. Avella and F. Mancini (Springer, Berlin, 2013) pp. 1-30. 
H. Lu, S. Sota, H. Matsueda, J. Bonca, and T. Tohyama, 
Phys. Rev. Lett. 109 , 197401 (2012). 

J. Rincon, K. A. Al-Hassanieh, A. E. Feiguin, and 

E. Dagotto, Phys. Rev. B 90 , 155112 (2014). 

L. B. Madsen, Phys. Rev. A 65, 053417 (2002). 

H. Castella, X. Zotos, and P. Prelovsek, Phys. Rev. Lett. 
74 , 972 (1995). 

X. Zotos and P. Prelovsek, Phys. Rev. B 53, 983 (1996). 

I. Bloch, J. Dalibard, and W. Zwerger, Rev. Mod. Phys. 
80 , 885 (2008). 

N. Tsuji, T. Oka, H. Aoki, and P. Werner, Phys. Rev. B 

85 , 155124 (2012). 


